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ABSTRACT 

A simple continuum damage mechanics (CDM) based 3D progressive damage 
analysis (PDA) tool for laminated composites was developed and implemented as a 
user defined material subroutine to link with a commercially available explicit finite 
element code. This PDA tool uses linear lamina properties from standard tests, 
predicts damage initiation with an easy-to-implement Hashin-Rotem failure criteria, 
and in the damage evolution phase, evaluates the degradation of material properties 
based on the crack band theory and traction-separation cohesive laws. It follows 
Matzenmiller et al.’s formulation to incorporate the degrading material properties 
into the damaged stiffness matrix. Since nonlinear shear and matrix stress-strain 
relations are not implemented, correction factors are used for slowing the reduction 
of the damaged shear stiffness terms to reflect the effect of these nonlinearities on the 
laminate strength predictions. This CDM based PDA tool is implemented as a user 
defined material (VUMAT) to link with the Abaqus/Explicit code. Strength 
predictions obtained, using this VUMAT, are correlated with test data for a set of 
notched specimens under tension and compression loads. 

INTRODUCTION 

Composites are lightweight and can be tailored to have superior stiffness and strength 
in loading directions. Thus, they have been widely used in primary aircraft structures. 
Composites are less damage tolerant than metals which can yield to redistribute loads. 
In addition, damage modes of composites are complex and can be difficult to detect. 
The building block approach has been used for development and certification of 
composite structures [1,2], however, this approach is time consuming and requires a 
large quantity of costly tests. To reduce the time and cost for development and 
certification of new composite aircraft structures, high fidelity progressive damage 
analysis (PDA) tools that can reliably predict the onset of damage and damage 
progression in composite structures are needed. 
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Many PDA tools have been developed and used [3-13]; however, accurately 
predicting the failure load of a composite structure is still a very challenging problem. 
A few traditional failure criteria have been widely used for predicting the damage 
initiation, such as the maximum stress criteria, the Hashin criteria [14], Hashin- 
Rotem criteria [15] and the Tsai-Wu [16] criteria. More recently developed failure 
criteria can be found in Refs. 17 to 19. There are various approaches for degrading 
ply properties in the damage evolution process [20]. For example, the stiffness of a 
damaged ply can be exponentially degraded, linearly degraded, or instantaneously 
degraded to zero [5-8]. The instantaneous degradation approach has been found to be 
too conservative and can predict lower structural strength than that determined by 
test data. Recent studies [9, 10] show that lamina property degradations after damage 
initiation can be modeled by using the crack band theory [9, 21]. The damage 
evolution process can thus be modeled as the opening of a cohesive crack. In the 
finite element analysis, the crack opening displacement is smeared into the damaged 
element using the characteristic length of the element resulting in a greater element 
strain, which in turn reduces the material property of the element. The characteristic 
length scale used in the crack band theory can also alleviate mesh dependency issues 
[22] often encountered in PDA. Furthermore, progression of damage in composite 
laminates often involves 3D stresses, including in-plane and transverse stresses. To 
include these stresses in the damage progression analysis, 3D PDA tools are required. 

In this study, a simple continuum damage mechanics (CDM) based 3D PDA tool is 
developed and implemented as a user defined material subroutine to link with a 
commercial finite element analysis software for laminated composite structures. The 
user defined material subroutine can be used as a framework for future development 
of more accurate and advanced PDA capabilities. The CDM approach [23-28] is well 
developed and easier to implement than discrete damage modeling [12, 13, 29, 30] 
approaches. In the CDM approach, damage is homogenized and treated as 
degradation of the material properties, thus no discrete cracks are modeled. This 
CDM based 3D PDA tool uses linear lamina properties from standard tests, predicts 
damage initiation with an easy-to-implement Hashin-Rotem failure criteria and in the 
damage evolution phase, evaluates the degradation of material properties based on 
the crack band theory and traction-separation cohesive laws. It follows Matzenmiller 
et al.’s formulation [28] to incorporate the degraded material properties into the 
damaged stiffness matrix. Since nonlinear shear and matrix stress-strain relations [9, 
18] are not implemented, correction factors are used for slowing the reduction of the 
damaged shear stiffnesses to reflect the effect of these nonlinearities on the laminate 
strength predictions. Furthermore, it implements length scales to alleviate the mesh 
dependency issue. A user defined material model is developed which links to a 
commercial explicit finite element analysis code for efficiently solving large linear 
and nonlinear composite structure problems. The solutions generated by this PDA 
tool are correlated with the experimental results. 
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This paper is organized as follows. First, the failure criteria used for determining the 
damage initiation in a lamina are presented. Second, the use of crack band theory and 
the traction-separation cohesive laws for determining the material property 
degradation in the damage evolution process are discussed and the conventional 
damage indices of all failure modes are defined. Third, a 3D CDM based model is 
established, using the damage indices as internal state variables. Fourth, the CDM 
equations for the development of user defined material model (VUMAT) for linking 
to Abaqus/Explicit [31] software are presented. Next, correlations of predictions with 
test data are presented. Last, concluding remarks are given at the end of the paper. 

FAILURE INITIATION CRITERIA 


For a general three-dimensional state of stress in terms of stress acting on the three 
principal material planes, the strain-based Hashin-Rotem failure criteria [15] used for 
this study are expressed as, 
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where s H and y.. ( i = 1,3 and j = 1,3 ) are the normal and engineering shear strains and 
the strain allowables A ( i = 1,6 ) are defined in Voigt notation as 
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The first failure criterion, Eq. 1, corresponds to the fiber failure, the second failure 
criterion, Eq. 2, corresponds to in-plane matrix dominated failure, and the third 
failure criterion, Eq. 3, corresponds to transverse matrix failure. 

DAMAGE EVOLUTION AND DEGRADATIONS OF MATERIAL 
PROPERTIES 


Once any of the Eqs. 1 to 3 is satisfied, the corresponding fiber and matrix damage 
initiates. The subsequent damage growth, namely, damage evolution, follows as the 
loading increases. The damage modes considered in this study are fiber tension and 
compression failure, matrix tension and compression failure, and the in-plane and 
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transverse shear failure. In this study, interlaminar delaminations are not modeled; 
the effects of delaminations on the strength predictions will be investigated in future 
studies. The evolution of each damage mode is assumed to follow a corresponding 
traction- separation law [9]. After damage initiates, it is assumed that the damage 
localizes and can be modeled as a cohesive crack within the element. As the element 
loading increases, the opening of the cohesive crack increases further. The crack 
opening displacement is smeared into the element, using the characteristic length of 
the element, resulting in increasing the element strain as shown in the Fig. 1. Using 
the crack band theory [21], the smeared strains of a damaged element are related to 
the crack opening displacements by [9] 
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where l e is a length scale for ply in-plane strains and l‘ e is the length scale for ply 
transverse strains. Note that sf, , s 22 , e 33 , yf 2 , yf 3 , and y 23 are the element damage 
initiation strains, determined by Eqs. 1 to 3, which remain the same for further load 
increments. Note that the shear strains yf 2 , yf 3 , and are engineering shear strains. 
In Eqs. 6-11, S { , 5 "' , S '" , S " 2 , S'" and <5™ are the crack opening displacements and 
any incremental change in the element strain after failure initiation is entirely from 
the crack opening displacement, where Sf , S 2 , S”‘ are normal, opening mode 
displacements, and S '" , S'” and S 23 are shearing mode displacements. Note that in 
this study, Abaqus 3D reduced integration element, C3D8R [31], is used, thus there 
is only one integration point for each element. Consequently, l e is the square root of 
the element in-plane area, and /' is the element thickness. The use of length scales 
can mitigate the dependence of the predictions on the mesh size [6, 9, 21, 22]. 

Pineda and Waas [9] derived the in-plane damaged moduli using the smeared strains 
and secant cohesive stiffnesses. The detailed derivations of the degraded 2D material 
properties can be found in Ref. 9. The derivations can be extended to predict degraded 
3D lamina properties for a laminate subject to monotonic and proportional loadings. 
The damaged tensile and shear moduli of a lamina can be expressed as 
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where E no , E 220 , E m , E m , E m , and E 230 are the linear elastic moduli of a 
unidirectional lamina and E n , E 22 , E 33 , G n , G n andG 23 are the degraded material 
moduli. G ft and G MT are the fiber and matrix tension fracture toughnesses, 
respectively, and G IIC is the Mode II critical energy release rate for shear failures. The 
t ( c , t 2C , t 3C , t"' 2c , f, 3C and t 23C are the maximum tractions of the triangular traction- 
separation laws [9] and they are assumed to be the stresses corresponding to the 
initiation strains determined by Eqs. 1 to 3. For material degradations under 
compression stresses, equations similar to Eqs. 12 to 17 are used. Note that the 
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unidirectional material properties for IM7/977-3 used in Eqs. 12 to 17 are given in 
Table I. 

Following the conventional definition [26, 28], these damage indices are defined as 
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Note that at an undamaged state (fully elastic), the index is set to be zero and at a 
fully damaged state, the index is set to be one. These damage indices are used as 
internal variables for establishing the 3D continuum damage mechanics model. 

3D CONTINUUM DAMAGE MECHANICS (CDM) MODEL 

Following the derivations in Ref. 28, the strain and true stress relation can be 
expressed as 
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The compliance matrix H can be expressed as, 
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Its inverse is the damaged stiffness matrix and the damaged stress-strain relation, i.e. 
the 3D CDM model, can be expressed as 
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C(d ) is the damaged stiffness matrix which can be expressed as 

C n C 12 C 13 0 0 0 

C 22 C 23 0 0 0 

C 33 0 0 0 

c 44 0 0 

Symmetry C 55 0 

C 6, 

where 

C 22 = (1 -d 2 )E 2 [1 - (1 - d x )(1 - d 3 ) v 13 v 31 ] / A 
C 33 =(l-d 3 )E 3 [l-(l-d 1 )(l-d 2 )v l2 v 21 ]/A 
C l2 — (1 — d t )(1 — d^)E ^\(\ — ^ 3 )^ 31^23 “t" v 21 ] / A 
C 13 = (1 - djXl - d 3 )E l [(1 - d 2 )v 21 v 32 + v 31 ] / A 
C 23 (1 d 2 ){ 1 <j? 3 )£' 2 [(1 )v 12 v 31 + v 32 ] / A 

C44 = (1 — d 4 )G 23 
C 55 = (1-c/ 5 )G 31 

C 66 =( 1 _ ^ 6) G 12 
and 

A = 1 — (1 — c/ 2 )(l — d 3 )\ 32 \ 23 -(l-dj)(l-d 3 )v 31 v 13 -(l-^ 1 )(l-^ 2 )v 21 v 12 
- 2(1 - dj)(l - d 2 )(l - d 3 )v 13 v 21 v 32 

Note that d t ( i = 1 , 6 ) are the damage variables which are functions of damage indices 
defined by Eqs. 18 to 23. The normal stiffness damaged variables di, d 2 , and di may 
be expressed as 
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Since the fiber damage, matrix damage, and shear damage modes can all affect the 
shear stiffness reduction, the shear damage variables d 4 , ds, and dr, may be expressed 
as the following phenomenological equations 

d 4 =l—(l-d 1 )(l-ad 2 )(l-ad 3 )(l-/3ds 23 ) (31) 

d 5 =1 -(1-^X1- ad 3 ){\- pds l3 ) (32) 

d 6 = 1 - (1-rf,)(1 - ad 2 )(l -/?ds 12 ) (33) 

Here, the values of a and ft are assumed to be between zero and one. When their 
values are zero, they do not affect the shear stiffness degradations. When their values 
are set to one, they have the maximum effect on the shear stiffness degradations. Note 
that the matrix and shear moduli for a unidirectional ply are assumed to be linear in 
this study. The a and ft may be regarded as correction factors to account for the 
effect of the nonlinearity of these moduli on the degradation of the damaged shear 
stiffnesses in Eq. 27. It is expected that nonlinear stress-strain curves for [90] n and 
[45/-45] ns can be used to calibrate the values of a and ft . Also, direct implementation 
of nonlinear stress-strain curves based on the Enhanced Schapery Theory (EST) [9] 
or the approach used in Ref. 1 8 may eliminate the need of using the correction factors 
a and p . 

USER DEFINED MATERIAL MODEL (VUMAT) FOR ABAQUS/EXPLICIT 

An ABAQUS/Explicit user defined material model (VUMAT) subroutine was 
developed to implement the CDM based 3D damaged constitutive equations. This 
VUMAT links to the Abaqus/Explicit code for PDA of laminated composite 
structures and is suitable for the reduced integration 3D element C3D8R [31]. 
Implementations for fully integrated 3D elements are possible, but will be more 
computationally intensive due to more integration points involved. In this VUMAT, 
the Green-Lagrange strain tensor is used, so it is applicable for material and 
geometrically nonlinear analyses. 

In the conventional CDM approach, the constitutive model is normally expressed in 
terms of stress-strain relations. When the material exhibits strain- softening behavior, 
leading to strain localization, this approach results in a strong mesh dependency of 
the finite element results due to the fact that energy dissipated decreases with mesh 
refinement [22]. In the crack band theory, the characteristic length scales are used, 
related to the element size, and the softening part of the constitutive law is expressed 
as a traction- separation relation. The use of characteristic length scales and traction- 
separation laws can assure that the energy released during the damage process is per 
unit area, not per unit volume, this can alleviate the mesh dependency issue [6, 9, 22]. 
The critical energy release rates (fracture toughnesses) are treated as additional 
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material parameters for determining the material property degradations in Eqs. 12 to 
17. In this VUMAT, the damage indices in Eqs. 18 to 23 are the solution-dependent 
state variables (SDVs) which each have a value between 0 and 1. A value of 0 means 
no damage while a value of 1 means total damage. These SDVs can be output for 
visualization to identify the extent of damage of each failure mode 

CORRELATIONS WITH EXPERIMENTAL DATA 

The prediction results for notched specimens were correlated with test data. The 
correlations can be used to verify if the VUMAT is properly coded and to evaluate 
the accuracy of the predictions. 

Experimental Data 

Experimental data of open hole tension (OHT) and open hole compression (OHC) 
specimens [11-13, 32-34] recently tested by the Air Force Research Laboratory 
(AFRL) were used for the test- analysis correlations. The dimensions of the specimens 
are shown in Fig. 2. For strain measurement, an extensometer was mounted on the 
specimen with its knife edges a half inch above and a half inch below the hole center, 
along the axial direction. These notched specimens were made of IM7/977-3 
laminates, containing various lay-up sequences, and were subjected to either tension 
or compression load [32, 33]. Three different layups were used for both OHT and 
OHC specimens; [0/45/90/-45]2s and [60/0/-60]3s were used as the strong layups 
which have 0°-plies and [30/60/90/-60/-30]2s was used as the weak layup which has 
no 0°-plies, to capture different failure modes. The cured average ply thickness was 
observed to be about 0.127 mm. These specimens were designed, fabricated and 
tested according to the ASTM standards. The unidirectional ply properties were also 
obtained using ASTM test standards. The linear elastic properties of a unidirectional 
IM7/977-3 ply are listed in Table I, which were obtained from Refs. 11-13 and 32- 
34. The fiber tension fracture toughness G FT and the fiber compression fracture 
toughness G FC shown in Table I were obtained from Refs. 7 and 35, respectively, for 
a similar graphite/epoxy system. The matrix tension fracture toughness G MT and the 
matrix compression fracture toughness G MC shown in Table I were obtained from 
Refs. 12, 13, and 34 for the IM7/977-3 material system. 

Finite Element Model 

The finite element model of the OHT and OHC specimens is shown in Fig. 3. 
Symmetric boundary conditions were used, and thus only the top half of the laminate 
was modeled. The number of nodes and the number of elements of the largest 
[30/60/90/-60/-30]2s model are 31,100, and 27,000, respectively. All the analyses 
were performed using high performance Linux computer clusters. The 
Abaqus/Explicit analyses were executed in double precision mode and the number of 
processors specified for each run was 12. The longest wall clock run time was about 
20 hours for the [30/60/90/-60/-30]2s OHT model and the shortest run time was about 
three hours for the [60/0/-60] 3 s OHC model. 
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Analysis Results 


In order to check the proper implementation of the VUMAT for elastic stiffness 
prediction of test specimens and to evaluate the sensitivity of the predicted test- 
specimen strength on the a and p values used, the stress and strain curves of the 
[45/0/90/-45]2 sOHT specimen with different a and p values are plotted in Fig. 4. A 
comparison with the test data shows that the elastic stiffness is well captured, thus 
verifying the implementation. In Fig. 4, the predictions using a- p = 0.3, 0.5 and 0.6 
are shown to illustrate that the changes of a and p values can affect the strength 
predictions. An increase of a and p values causes greater shear stiffness reductions 
in the damage evolution stage which in turn can lower the predicted strengths. Note 
that the change of a and p values from 0.3 to 0.6 reduces the predicted strength by 
only 7%, indicating that the predicted strength is not be not very sensitive to the a 
and p values used. Similar trends were observed for other layups. For the following 
tested and predicted strength correlations performed for the OHT and OHC 
specimens, all the analysis results were obtained with a - P-0.5. It was observed 
that this set of a and p values served as suitable correction factors to slow the shear 
stiffness reductions for all the cases analyzed and peak loads were reached without 
encountering any large element distortion problem. 

The predicted strengths of the OHT and OHC specimens with three different lay-ups 
are compared with the experimental data as shown in Table II. The prediction errors 
with the averaged test data are also presented in the table. These strengths were 
obtained by dividing the peak load with the total cross-section area. The total cross- 
section area is the specimen width multiplied by the laminate thickness. For example, 
the peak load for the quasi-isotropic OHT specimen can be clearly identified from 
load displacement curve shown in Fig. 5. The predicted strengths and test data are 
also plotted in Fig. 6 for visual comparisons. The analysis results presented in Table 
II and Fig. 6 correlate reasonably well with test results, considering the simplicity of 
this VUMAT implementation. However, some errors are still more than 10%, 
indicating further development or adoption of advanced PDA tools are needed for 
more accurate predictions. Since the interlaminar delamination was not modeled in 
this study, the effect of delamination on the strength prediction needs be investigated 
in future studies. 

To illustrate that all specimens analyzed using a = p = 0.5 do reach their peak loads, 
indicated by the accumulation of significant damage across the width of the 
specimens, the fiber tensile failure for the outermost major load bearing ply (see red 
arrow) for each OHT specimen is shown in Fig. 7 and the fiber compression failure 
for the outermost major load bearing ply for each OTC specimen is shown in Fig. 8. 
These figures indicate that the major load bearing plies are nearly failed after the peak 
load has been reached. Note that a complete fiber tensile failure is indicated by SDV 1 
reaching a value of 1.0 in Fig. 7 while a complete fiber compression failure is 
indicated by SDV2 reaching a value of 1.0 in Fig. 8. 
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CONCLUDING REMARKS 


A simple CDM based 3D PDA tool has been developed and implemented as a user 
defined material (VUMAT) subroutine to link with a commercial finite element code, 
ABAQUS/Explicit. This PDA tool uses linear lamina properties from ASTM 
standard tests, predicts damage initiation with an easy-to-implement Hashin-Rotem 
failure criteria, and in the damage evolution phase, evaluates the degradation of 
material properties based on the crack band theory and traction- separation cohesive 
laws. It follows Matzenmiller et al.’s formulation to incorporate the degrading 
material properties into the damaged stiffness matrix. The user defined material 
subroutine can serve as a framework for future development of more accurate and 
advanced PDA capabilities. 

The characteristic length scales related to the crack band theory were used to alleviate 
the mesh dependency issue often associated with PDA predictions. This VUMAT 
uses phenomenological parameters to control the shear stiffness degradations. Note 
that the matrix and shear moduli for a unidirectional ply are assumed to be linear in 
this study. The a and (j may be regarded as correction factors to account for the 
effect of the nonlinearity of these moduli on the degradation of the damaged shear 
stiffnesses. In future studies, the a and (j values may be judiciously calibrated by 
using the nonlinear matrix and shear stress-strain curves obtained from simple 
coupon tests. Preliminary predictions of the strengths were in reasonably good 
agreement with test data of OHT and OHC specimens with various layups when 
a = p = 0.5 was used. However, some errors are still more than 10%, indicating the 
need for modeling interlaminar delaminations and further development or adoption 
of advanced PDA tools for more accurate predictions. 
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Table I - MATERIAL PROPERTIES 


Material Property 


E no (GPa ) 

164.0 (Tension) 
137.4 (Compression) 

E 220 5 E;, 1(l ( GPa ) 

8.98 

G m’G m (GPa) 

5.02 

G 230 (GPa ) 

3.00 

Vl2.V 13 

0.320 

V 23 

0.496 

G IC ,G IIC (N / mm) 

0.256, 1.156 

Y T ,Y c (MPa) 

100.0, 247.0 

S(MPa) 

80.0 

X T ,X c {GPa) 

2.90, 1.68 

G FT ,G FC (N / mm) 

81.534, 24.533 

G MT ’ ^MC ^ tWfl) 

0.256, 1.156 


Table II - TEST AND ANALYSIS RESULTS 


Lay-Ups 

OHT 

OHC 

Analysis 

(MPa) 

Test 

(MPa) 

% Error 

Analysis 

(MPa) 

Test 

(MPa) 

% Error 

[0/45/90/-45l2s 

508 

554 

-9.06 

318 

341 

-6.74 

f60/0/-601ss 

533 

543 

-1.88 

290 

358 

-18.99 

[30/60/90/-60/-30l2s 

461 

409 

11.28 

327 

295 

10.85 
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Figure 1. Damaged element with cohesive opening. 


138 mm ► 



Clamped 

Edge 


Figure 2. Dimensions of IM7/977-3 open hole specimen. 
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Figure 3. 3D finite element mesh of open hole specimen. 



Figure 4. Stress-strain plot for OHT [0/45/90/-45]2s specimen using various 

values of a and p . 


16 



Stress (MPa) % Force < N ) 



5. Predicted load-displacement curve for [0/45/90/-45]2s OHT specimen. 



OHT OHC OHT OHC OHT OHC 
[0/45/90/-45] 2S [60/0/-60] 3S [30/60/90/-60/-30] 2S 

Figure 6. Test and analysis results correlations. 
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Figure 7. Final fiber tensile failure (SDVI) in outermost major load bearing ply. 
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Figure 8. Final fiber compression failure (SDV2) in outermost major load bearing 

ply- 


is 







